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LONG-TERM  GOALS 

The  long-term  goal  of  this  project  is  to  design  and  evaluate  the  components  that  will  comprise  a  next 
generation  mesoscale  atmospheric  model  within  the  Coupled  Ocean/ Atmosphere  Mesoscale  Prediction 
System  (COAMPS®1).  It  is  anticipated  that  in  order  to  meet  future  Navy  requirements,  next  generation 
approaches  to  numerical  techniques  and  physical  parameterizations  will  be  needed. 

OBJECTIVES 

The  objectives  of  this  project  involve  the  development,  testing,  and  validation  of  new  numerical  tech¬ 
niques  such  as  spectral  element  techniques  and  monotonic  advection  methods.  The  overarching  objec¬ 
tive  is  to  develop  methods  for  high-resolution  numerical  weather  prediction  applications  for  horizontal 
grid  increments  at  1  km  or  less. 

APPROACH 

Our  approach  is  to  follow  a  methodical  plan  in  the  development  and  testing  of  a  nonhydrostatic  micro¬ 
scale  modeling  system  that  will  leverage  the  existing  COAMPS  and  new  model  prototypes.  Our  work 
on  numerical  methods  will  involve  investigation  of  spatial  and  temporal  discretization  algorithms  that 
are  superior  to  the  current  generation  leap-frog,  second-order  accurate  numerical  techniques  presently 
employed  in  COAMPS  and  many  other  models;  these  new  discretization  methods  will  be  developed 
and  implemented.  We  would  like  to  use  numerical  methods  that  have  already  been  developed  and 
tested  in  other  communities,  such  as  computational  fluid  dynamics,  and  apply  these  to  high-resolution 
numerical  weather  prediction  applications.  Validation  and  evaluation  of  the  modeling  system  will  be 
performed  using  datasets  of  opportunity,  particularly  in  regions  of  Navy  significance. 


1  CO  AMPS 8  is  a  registered  trademark  of  the  Naval  Research  Laboratory. 
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WORK  COMPLETED 


1.  Spectral  element  2D  prototypes. 

Spatial  resolution  required  to  adequately  resolve  features  of  interest  can  be  easily  determined  for  the 
finite  difference  model  where  the  grid  spacing  is  constant.  Even  with  varying  resolution  (e.g.  common¬ 
ly  used  vertically  stretched  grid)  the  grid  spacing  is  increasing  monotonically  with  the  distance  away 
from  the  area  of  interest.  The  non-uniform  spacing  of  nodal  points  used  for  the  spectral  element  model 
makes  the  same  task  less  straightforward.  The  difference  between  the  minimum  and  maximum  dis¬ 
tance  between  adjacent  nodal  points  within  the  same  element  is  proportional  to  the  polynomial  order. 
The  h-p  parameter  space,  where  h  is  number  of  elements  in  the  horizontal  direction  and  p  is  polynomi¬ 
al  order  was  mapped  out  to  span  nominal  horizontal  resolution  from  200  to  10000  m. 

Numerical  solutions  to  the  linear,  hydrostatic  flow  over  an  idealized  mountain  and  an  analytic  solution 
for  the  vertical  momentum  flux  were  used  to  calculate  the  normalized  1 2  norm.  An  overall  error,  speed 
of  convergence  and  computational  cost  were  assessed,  the  latter  two  compared  to  the  numerical  solu¬ 
tions  obtained  with  a  finite  difference  model. 

There  is  no  analytical  solution  for  our  chosen  test  case  involving  moisture  (squall  line  with  warm  mi¬ 
crophysics).  Storm  stages,  checked  at  two  different  times  (3000  and  6000  s)  differ  depending  on  the 
polynomial  order  and  nominal  grid  spacing.  In  addition,  dimensional  viscosity  was  added  to  the  model 
to  ensure  stability  and  realism  of  the  solution  since  there  is  no  subgrid-scale  mixing.  The  initial  value 
of  viscosity  was  constant  in  all  cases,  regardless  of  the  time  step  chosen.  Tests  with  various  values  of 
viscosity  were  perfonned  to  assess  convergence  of  solutions.  Symmetry  tests  were  perfonned  to  assess 
sensitivity  of  the  storm  development  to  the  initial  warm  bubble  location  in  relation  to  the  non- 
uniformly  spaced  nodal  points. 

2.  WENO  (Weighted  Essentially  Non-Oscillatory)  methods 

Atmospheric  models  require  numerical  methods  that  can  accurately  represent  the  transport  of  tracers 
with  steep  gradients,  such  as  those  that  occur  at  cloud  boundaries  or  the  edges  of  chemical  plumes.  In 
atmospheric  sciences,  the  most  widely  used  numerical  techniques  for  this  type  of  problem  are  flux- 
corrected  transport  or  closely  related  flux-limiter  methods.  The  limiters  are  typically  designed  to  pre¬ 
vent  the  development  of  new  extrema  in  the  concentration  field.  This  will  preserve  the  non-negativity 
of  initially  non-negative  fields,  which  is  essential  for  the  correct  simulation  of  cloud  microphysics  or 
chemical  reactions.  One  serious  systematic  weakness  of  flux  limiter  methods  is  that  they  also  tend  to 
damp  the  amplitude  of  extrema  in  smooth  regions  of  the  flow,  such  as  the  trough  of  a  well-resolved 
sine  wave.  To  avoid  this  problem,  we  have  been  investigating  the  application  of  WENO  (Weighted 
Essentially  Non-Oscillatory)  methods  to  tracer  transport  in  atmospheric  models.  WENO  methods  are 
widely  used  in  many  disciplines,  but  scarcely  been  tested  in  atmospheric  applications.  WENO  me¬ 
thods  preserve  steep  gradients  while  simultaneously  avoiding  the  dissipation  of  smooth  extrema  by  es¬ 
timating  the  value  of  the  solution  in  a  way  that  heavily  weights  the  smoothest  possible  cubic  poly¬ 
nomial  fit  to  the  local  function  values.  Where  the  solution  is  well  resolved,  all  possible  cubic  interpo- 
lants  are  weighted  almost  equally.  Near  a  steep  gradient,  those  interpolants  are  almost  completely  ig¬ 
nored. 
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RESULTS 


1.  Spectral  element  2D  prototypes. 

The  normalized  12  norm  error  is  lowest  for  the  nominal  resolution  less  than  or  equal  to  1  km  (Figure 
1).  Large  error  in  the  lower  left  portion  of  the  same  figure  is  due  to  the  poorly  resolved  topography  and 
error  introduced  with  inexact  integration  when  the  polynomial  order  is  too  low.  The  convergence  to  the 
final  solution  is  achieved  quickly  for  cases  with  coarser  nominal  resolution  of  2000  and  3000  m  (yel¬ 
low  and  green  lines,  respectively  in  Figure  2),  while  simulations  with  finer  nominal  resolution  con¬ 
verge  at  or  around  12  hours.  The  most  important  factor  in  the  final  error  is  the  nominal  horizontal  grid 
spacing,  where  all  the  cases  can  be  grouped  into  two  clusters:  quickly  converging  with  less  accuracy 
and  slowly  with  more  accuracy.  Note  that  the  errors  calculated  with  the  finite  difference  model  keep 
decreasing  with  increasing  resolution,  but  even  with  the  finest  resolution  (500  m,  solid  gray  line  with 
diamonds)  the  error  is  comparable  to  the  error  obtained  with  the  spectral  element  model  with  the  no¬ 
minal  resolution  of  2000  m. 

The  finite  difference  model  is  computationally  cheaper  for  a  given  nominal  resolution  (Figure  3).  The 
increase  of  computational  cost  per  increased  resolution  is  comparable  for  both  models.  While  the  reso¬ 
lution  refinement  results  in  a  monotonic  decrease  of  the  error  for  the  finite  difference  model,  thus  justi¬ 
fying  the  cost,  which  is  not  the  case  for  the  spectral  element  model.  Increasing  the  nominal  resolution 
from  3000  to  2000  m,  or  1000  to  500  m,  yields  only  marginal  error  reduction,  with  an  increased  com¬ 
putational  cost.  The  biggest  accuracy  yield  happens  when  reducing  the  average  resolution  from  2000 
to  1000  m  (more  than  an  order  of  magnitude). 

The  nominal  horizontal  resolution  of  2000  m  is  too  coarse,  or  the  effective  diffusion  combined  with  a 
larger  time  step  is  too  weak,  resulting  in  larger  storm  clouds  with  weaker  positive  virtual  potential 
temperature  perturbations  compared  to  those  obtained  with  finer  resolutions  (Figure  4,  right  panel). 
More  consistent  results  were  simulated  with  both  500  and  1000  m  nominal  resolutions  (Figure  4,  left 
and  middle  panel,  respectively). 

Smaller  time  step  required  by  finer  resolution  results  in  applying  diffusion  more  often.  In  order  to  bal¬ 
ance  the  effective  diffusion,  it  was  increased  for  1000  m  (larger  time  step  compared  to  500  m)  compar¬ 
ing  the  results  (Figure  5,  right  panel)  to  the  original  500  m  (Figure  4,  left  panel).  The  narrower  updraft 
section  of  the  cloud  is  now  more  comparable.  Similarly,  a  decrease  in  diffusion  for  the  500  m  simula¬ 
tion  resulted  in  a  broader  updraft  section  (Figure  5,  left  panel),  compared  to  the  original  1000  m  case 
(Figure  4,  middle  panel). 

The  initial  location  of  the  wann  bubble,  which  triggers  the  storm,  is  important  for  the  symmetry  of  the 
fully  developed  storm.  When  the  bubble  is  centered,  such  that  the  nodal  spacing  to  the  left  and  right 
are  equal,  the  storm  will  remain  symmetric  in  the  absence  of  wind  shear.  If  the  nodal  spacing  is  not 
equal,  the  storm  develops  a  tilt,  similar  to  simulations  with  the  environmental  wind  shear  (not  shown). 

2.  WENO  (Weighted  Essentially  Non-Oscillatory)  methods 

The  WENO  based  selective  monotonic  advection  (SMA)  scheme  (Blossey  and  Durran  2007)  was  fully 
implemented  in  COAMPS  in  FY  2009.  In  addition,  a  semi-Lagrangian  option  (Skamarock  2006;  Blos¬ 
sey  and  Durran,  2007)  was  also  implemented  in  order  to  offset  the  extra  computational  costs  required 
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for  the  SMA  scheme.  Both  the  semi-Lagrangian  and  Eulerian  versions  of  the  SMA  scheme  were 
tested  and  compared  with  the  currently  available  2nd-  and  4th-order  finite  difference  advection  schemes 
for  a  variety  of  idealized  test  problems.  In  addition,  both  SMA  schemes  were  evaluated  in  the  full 
COAMPS  numerical  weather  prediction  (NWP)  mode  over  two  15-day  periods:  1-15  January  2008  for 
an  eastern  Pacific  domain  and  1-15  May  2008  for  a  domain  containing  the  continental  United  States 
(CONUS).  The  CONUS  test  was  configured  with  a  9-km  horizontal  resolution  domain  extending  from 
the  Rocky  Mountains  to  the  East  Coast  nested  within  a  27-km  horizontal  resolution  domain  containing 
the  entire  CONUS  region.  Several  highlights  from  the  9-km  domain  will  be  discussed  below. 

One  issue  with  any  finite  difference  scheme  is  the  presence  of  overshoots  and  undershoots  in  regions 
of  steep  gradients.  For  fields  that  are  required  to  be  positive  definite,  such  as  moisture,  any  negative 
values  generated  must  be  filled  in  order  to  maintain  a  physically  meaningful  solution.  Skamarock  and 
Weisman,  2006  have  shown  that  the  cumulative  effect  of  filling  negative  values  of  moisture  can  lead  to 
a  positive  bias  in  precipitation.  They  demonstrate  that  a  positive  definite  advection  improves  this  bias 
by  eliminating  the  necessity  to  fill  negative  values.  The  impact  of  the  positive  definite  SMA  scheme 
on  the  COAMPS  precipitation  bias  is  shown  in  Fig.  6.  Plotted  is  the  24-h  accumulated  precipitation 
model  bias  as  a  function  of  the  accumulated  precipitation  for  the  fourth-order  advection  (black),  Eule¬ 
rian  SMA  scheme  (blue),  and  semi-lagrangian  SMA  scheme  (red).  Here  the  time  step  of  the  semi- 
lagrangian  scheme  is  three  times  that  of  the  Eulerian  scheme.  The  24-  and  48-  hr  forecasts,  sampled 
from  the  9-km  CONUS  domain  for  the  1-15  May  2008  time  period,  are  compared  with  the  Stage-IV 
gridded  precipitation  data  set  (Lin  and  Mitchell  2005).  Bias  values  less  than  1  indicate  a  negative  bias, 
1  is  no  bias,  and  values  greater  than  1  indicate  a  positive  bias.  While  the  each  scheme  shows  and  in¬ 
creasing  bias  as  a  function  of  accumulated  precipitation,  both  the  Eulerian  and  semi-Lagrangian  SMA 
schemes  have  a  lower  bias  than  the  fourth-order  advection  scheme.  This  is  especially  true  for  larger 
precipitation  thresholds. 

Figure  7  compares  the  morphology  of  the  2-km  AGL  cloud- water  mixing  ratio  field  for  a  pair  of 
COAMPS  simulations  taken  from  the  9-km  CONUS  domain  and  initialized  00  UTC  14th  May  2008. 
This  period  was  chosen  due  to  the  presence  of  a  developing  cyclone,  trailing  cold  front,  and  active 
convection  over  the  Midwest.  As  the  cold  front  progressed  through  Texas,  Oklahoma,  and  Arkansas 
and  into  the  Gulf  Coast  states,  radar  imagery  indicated  several  developing  lines  of  severe  thunders¬ 
torms  (not  shown).  Both  the  fourth-order  advection  scheme  and  Eulerian  SMA  scheme  indicate  this 
line  of  thunderstorms  at  00  UTC  15th  May  (24-hr  forecast;  Fig.  7a  and  b).  However,  over  the  next  12 
hours,  as  the  line  progresses  into  Louisiana  and  Mississippi,  the  SMA  scheme  is  better  able  to  maintain 
a  coherent  linear  structure  in  the  cloud  field  (Fig.  7d)  compared  to  the  fourth-order  advection  scheme 
were  the  cloud  field  is  much  more  disjoint  and  consists  of  individual  blobs  of  cloud-liquid  water  (Fig. 
7c).  The  cloud  field  present  in  the  simulation  with  the  SMA  scheme  is  more  consistent  with  the  preci¬ 
pitation  field  from  the  observed  radar  imagery  (not  shown);  however,  an  objective  verification  compar¬ 
ison  of  the  cloud  fields  in  the  two  schemes  has  not  been  completed. 

IMPACT/APPLICATIONS 

COAMPS  is  the  Navy’s  operational  mesoscale  NWP  system  and  is  recognized  as  the  key  model  com¬ 
ponent  driving  a  variety  of  DoD  tactical  decision  aids.  Accurate  mesoscale  prediction  is  considered  an 
indispensable  capability  for  defense  and  civilian  applications.  Skillful  COAMPS  predictions  at  resolu¬ 
tions  less  than  1  km  will  establish  new  capabilities  for  the  support  of  the  warfighter  and  Sea  Power  2 1 . 
Operational  difficulties  with  weapon  systems  such  as  the  Joint  Standoff  Weapon  (JSOW)  have  been 
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documented  in  regions  with  fine-scale  topography  due  to  low-level  wind  shear  and  turbulence.  Im¬ 
proved  high-resolution  predictive  capabilities  will  help  to  mitigate  these  problems  and  introduce  poten¬ 
tially  significant  cost  saving  measures  for  the  operational  application  of  JSOW.  The  capability  to  pre¬ 
dict  the  atmosphere  at  very  high  resolution  will  further  the  Navy  sea  strike  and  sea  shield  operations, 
provide  improved  representation  of  aerosol  transport,  and  will  lead  to  tactical  model  improvements. 
Emergency  response  capabilities  and  Homeland  Security  issues  within  the  DoD  and  elsewhere,  such  as 
LLNL,  will  be  enhanced  with  the  new  modeling  capability. 

TRANSITIONS 

The  next  generation  COAMPS  system  will  transition  to  6.4  projects  within  PE  0603207N  (SPAWAR, 
PMW-120)  that  focus  on  the  transition  COAMPS  to  FNMOC.  The  improvements  to  the  COAMPS 
dynamical  core  have  been  transitioned  to  the  SPAWAR  6.4  project  and  subsequently  to  operations  as  a 
result  of  the  marked  improvement  in  the  geopotential  height  bias  statistics. 

RELATED  PROJECTS 

COAMPS  will  be  used  in  related  6. 1  projects  within  PE  0601 153N  that  include  studies  of  air-ocean 
coupling,  boundary  layer  studies,  and  topographic  flows  and  in  related  6.2  projects  within  PE 
060243  5N  that  focus  on  the  development  of  the  atmospheric  components  (QC,  analysis,  initialization, 
and  forecast  model)  of  COAMPS. 

REFERENCES 

Blossey,  P.N.  and  D.R.  Durran,  2007:  A  simple,  effective  WENO-like  smoothness  metric  for  use  in 
conservative  models  of  scalar  advection.  J.  Comput.  Phys.,  227,  5160-5183. 

Lin,  Y.  and  K.  E.  Mitchell,  2005:  The  NCEP  stage  II/IV  hourly  precipitation  analyses:  Development 
and  applications.  Preprints,  19th  Conf.  on  Hydrology,  San  Diego,  CA,  Amer.  Meteor.  Soc.,  CD- 
ROM,  1.2. 

Skamarock,  W,  2006:  Positive-definite  and  monotonic  limiters  for  unrestricted-time-step  transport 
schemes.  Mon.  Wea.  Rev.,  134,  2241-2250. 

PUBLICATIONS 

Doyle,  J.D.,  Q.  Jiang,  R.B.  Smith,  V.  Grubisic,  2010:  Three-dimensional  characteristics  of  stratos¬ 
pheric  mountain  waves  during  T-REX.  Mon.  Wea.  Rev.  (In  Press) 

Doyle,  J.D.,  C.C.  Epifanio,  A.  Persson,  P.A.  Reinecke,  G.  Zangl,  2010:  Mesoscale  modeling  over 
complex  terrain:  Numerical  and  predictability  perspectives.  AMS  Monograph  Complex  Terrain 
Processes  and  Forecasting.  (In  Press) 

Jiang,  Q.,  M.  Liu,  and  J.  D.  Doyle,  2010:  Influence  of  Mesoscale  Dynamics  and  Turbulence  on  Fine 
Dust  Transport  in  Owens  Valley.  J.  Applied  Meteor.  Clim.  (In  Press) 

Schmidli,  J.,  B.  J.  Billings,  F.  K.  Chow,  J.  D.  Doyle,  V.  Grubisic,  T.  R.  Holt,  Q.  Jiang,  K.  A.  Lund- 
quist,  P.  Sheridan,  S.  Vosper,  S.  F.  J.  De  Wekker,  C.  D.  Whiteman,  A.  A.  Wyszogrodzki,  G. 

Zaengl,  2011:  Intercomparison  of  mesoscale  model  simulations  of  the  daytime  valley  wind  sys¬ 
tem,  Mon.  Wea.  Rev.,  In  Press. 


5 


0.8 

0.7 

0.6 

0.5 

0.4 

0.3 

0.2 

0.1 

6  8  10  12  16  20  24  30  40  48  60  80  120 

/^refinement  (number  of  elements  in  horizontal) 


o> 

T3 


E 

o 

C 

o 

Q. 


c 

E 

<D 

C 

V 

s_ 

in 


Figure  1.  Normalized  1 2  norm  distribution  (shaded  contours,  c.i.  0.1)  of  the  vertical  momentum 
flux  as  a  function  of  the  polynomial  order  (p)  and  number  of  elements  in  horizontal  (h)  . 
Red,  orange,  yellow  and  green  lines  connect  cases  with  the  same  average  resolution  of 500, 

1000,  2000  and  3000  m,  respectively. 
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Figure  2.  Time  evolution  of  the  normalized  1 2  norm  of  the  vertical  momentum  flux  (output  every 
hour,  starting  at  t=l  h).  Results  based  on  simulations  with  the  same  horizontal  resolution  are 
grouped  by  a  line  color:  green  (3.0  km),  yellow  (2.0  km),  orange  (1.0  km)  and  red  (0.5  km).  Line 
styles  depict  the  polynomial  order  of  basis  functions:  thick  solid  (p=4),  dotted  (p=6),  dashed  (p=8) 
and  thin  solid  (p=10).  In  addition,  gray  lines  with  diamonds  represent  results  obtained  with  a  finite 
difference  model:  solid,  short  dashed,  long  dashed  and  dotted  stand  for  grid  spacing  of 500,  1000, 

2000  and  3000  m,  respectively. 
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Figure  3.  Normalized  1 2  norm  of  the  vertical  momentum  flux  as  a  function  of  computational  time. 
Results  obtained  with  the  spectral  element  model  are:  solid,  dashed,  dot-dashed,  and  dotted  lines  for 
p=4,  p=6,  p=8  and  p=10,  respectively.  The  lightest  blue  line  with  triangles  represents  the  finite  dif¬ 
ference  model.  Simulations  with  different  resolutions  are  represented  with  small  red  circles  (500  m), 
orange  circles  with  wide  rings  (1000  m),  yellow  circles  with  thick  inner  and  thin  outer  rings  (2000 

m)  and  green  circles  with  two  thick  rings  ( 3000  m). 


Figure  4.  Squall  line  evolution  at  t=6000  s  for  the  nominal  resolution  of 500,  1000  and  2000  m, 
in  left,  middle  and  right  panel,  respectively.  The  cloud  outline  is  based  on  the  mixing  ratio  of 
cloud  water  (10~5),  colored  contours  are  virtual  potential  temperature  perturbations  (c.i.  3  K, 
negative/positive  values  are  blue/red).  In  addition,  storm  relative  circulation  stronger  than  1  m/s  is 
depicted  by  arrows,  blue/red  for  sinking/rising.  Dimensional  diffusivity  in  all  cases  was  200  m2/s, 
time  step  was  0.2,  0.25  and  0.5  s  for  cases  with  nominal  horizontal  resolution  of 500,  1000  and 

2000  m,  respectively. 
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Figure  5.  Same  as  Figure  4,  except  with  weaker  diffusivity  (1 60  m2/s)  for  the  case  with  nominal  res¬ 
olution  of 500  m  (left  panel)  and  stronger  diffusivity  (250  m/s)  for  1000  m  (right  panel). 


24-h  Forecast  Bias 


Threshold  (mm) 


48-h  Forecast  Bias 


Figure  6.  The  COAMPS  model  bias  of  24-h  accumulated  precipitation  from  24-  and  48-h  forecasts 
during  May  2008.  The  horizontal  resolution  is  9-km  and  the  computational  covers  the  continental 

United  States  east  of  the  Rocky  Mountains. 
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Figure  7.  A  pair  COAMPS  model  forecasts  for  the  (a,  c)  4th -order  advection  scheme  and  the  (b,  d) 
WENO-based  selective  monotonic  advection  scheme.  Plotted  is  the  (a-b)  24-  and  (c-d)  36-h  forecast 
of  the  2-km  cloud-water  mixing  ratio  (g/kg)  in  a  subset  of  a  9-km  horizontal  resolution  domain.  The 

model  forecast  was  initialized  00  UTC,  14  May. 
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